%First plotting for panels 1B, 1C

load('dataABR')
imagingAnimalsNE={'AEH147';
    'AEH155';
    'AEH162';
    'AEH165';
    'AEH168';
    'AEH170';
    'AEH174';
    'AEH178'};

HLexploreAnimals={'AEH141';
    'AEH144';
    'AEH145';
    'AEH149';
    'AEH150'};

behaviorAnimalsNE={'MM06';
    'MM07';
    'MM09';
    'MM12';
    'MM18';
    'MM20';
    'MM22';};

frequencyValues=[4 5.7 8 11.3 16 22.6 32 45.3 64]; 
imgDays={'day_pre';'day_1';'day_3';'day_5';'day_7'}; nDays=size(imgDays,1); colorDays=copper(nDays);
[ABRthr,ABRstdE,ABRn]=concatenateABRs(dataABR,[HLexploreAnimals;behaviorAnimalsNE;imagingAnimalsNE],frequencyValues,imgDays);

figure
for dd=1:nDays
    shadedErrorBar(frequencyValues.*1000,squeeze(ABRthr(:,dd)),squeeze(ABRstdE(:,dd)),'lineprops',{'color',colorDays(dd,:),'linewidth',3})
    hold on
    p{dd} = plot(frequencyValues.*1000,squeeze(ABRthr(:,dd)),'color',colorDays(dd,:),'linewidth',3);
end
axis square
box off
ylim([10 105]);
ylabel('ABR wave1 threshold (dB SPL)')
xlabel('Tone Frequency (kHz)')
set(gca,'fontsize',14)
set(gca,'XScale','log')
xlim([4000 64000])
SetLogTicks('x',2,2000,64000)
legend([p{1} p{2} p{3} p{4} p{5}],{'Baseline' 'Day 1' 'Day 3' 'Day 5' 'Day 7+'})



figure
imgDays={'day_pre';'day_1';'day_3';'day_5';'day_7'}; nDays=size(imgDays,1); colorDays=copper(nDays);
[DPOAEthr,DPOAEstdE,DPOAEn]=concatenateDPOAEs(dataABR,[HLexploreAnimals;behaviorAnimalsNE;imagingAnimalsNE],frequencyValues,imgDays);
hold off
for dd=1:nDays
    shadedErrorBar(frequencyValues.*1000,squeeze(DPOAEthr(:,dd)),squeeze(DPOAEstdE(:,dd)),'lineprops',{'color',colorDays(dd,:),'linewidth',3})
    hold on
    p{dd} = plot(frequencyValues.*1000,squeeze(DPOAEthr(:,dd)),'color',colorDays(dd,:),'linewidth',3);
end
axis square
box off
ylim([10 105]);
ylabel('DPOAE threshold (dB SPL)')
xlabel('Tone Frequency (kHz)')
set(gca,'fontsize',14)
set(gca,'XScale','log')
xlim([4000 64000])
SetLogTicks('x',2,2000,64000)
legend([p{1} p{2} p{3} p{4} p{5}],{'Baseline' 'Day 1' 'Day 3' 'Day 5' 'Day 7+'})

%%
%Statistics for Figure 1B,1C Let's start by first looking at the ABR's. We'll
%do this with a repeated measures ANOVA, with the effect of frequency and
%of time (but here only using two time points: baseline and 7+ post
%exposure time point).

%First set the names for the mice
clear; clc;
load('dataABR')
imagingAnimalsNE={'AEH147';'AEH155';'AEH162';'AEH165';'AEH168';'AEH170';'AEH174';'AEH178'};
HLexploreAnimals={'AEH141';'AEH144';'AEH145';'AEH149';'AEH150'};

%Now set the values
frequencyValues=[4 5.7 8 11.3 16 22.6 32 45.3 64]; 
imgDays={'day_pre';'day_1';'day_3';'day_5';'day_7';'day_10'}; 
nDays=size(imgDays,1); animals = [HLexploreAnimals;imagingAnimalsNE];
nFreq=length(frequencyValues); nDays=length(imgDays); nAnimals=length(animals);

%Store in a 3D matrix that is animals x frequency x day
thrABR=nan(nAnimals,nFreq,nDays);

%Create the matrix
for ani=1:nAnimals
    animalStr=animals{ani,:};
    for dd=1:nDays
        dayStr=imgDays{dd,:};
        if ~isfield(dataABR.(animalStr),dayStr)
            continue
        end
        for freq=1:nFreq
            freqStr=['freq_' num2str(round(frequencyValues(freq)))];
            if isfield(dataABR.(animalStr).(dayStr).both,freqStr)
                thrABR(ani,freq,dd)=dataABR.(animalStr).(dayStr).both.(freqStr).threshold;
            end
        end
    end
end

%Now we're going to select two specific time points, and only use
%frequencies that have been tested across all mice.

for i = 1:size(thrABR,1) %Loop through mice
    nums1 = thrABR(i,1:9,1); %Get baseline measurements
    nums2 = nan(1,9); day = size(thrABR,3);
    while any(isnan(nums2(3:8)))
        nums2 = thrABR(i,1:9,day);
        day = day - 1;
    end
    allValues(i,:) = [nums1 nums2];
end

%Now we want to format this into a table with the correct labels
for i = 1:18 %9 frequencies, 2 time points
    v = strcat('V',num2str(i));
    varNames{i,1} = v;
end

dataTable = array2table(allValues, 'VariableNames', varNames);

%Finally create a table that reflects the within subject factors, which
%here are frequency and time. 
factorNames = {'Frequency'; 'Day'};

freqLabels = repmat(frequencyValues,1,2)'; 
freqLabels = arrayfun(@num2str, freqLabels, 'UniformOutput', 0);

dayLabels = [ones(1,9) ones(1,9).*2]';
dayLabels = arrayfun(@num2str, dayLabels, 'UniformOutput', 0);

withinTable = table(freqLabels,dayLabels,'VariableNames',factorNames);

%Now that we have the tables set up, run the anova model
rmABR = fitrm(dataTable, 'V1-V18~1','WithinDesign',withinTable);
[rAnovaResults] = ranova(rmABR, 'WithinModel','Frequency*Day');
rAnovaResults

%Also do post-hoc tests for each frequency. Do as a paired t test, times
%two to correct for multiple comparisons.
for i = 1:9 %For 9 frequencies
    nums1 = allValues(:,i); nums2 = allValues(:,(i+9));
    [h,p] = ttest(nums1,nums2);
    disp(freqLabels{i}); disp(p*2); disp('');
end

%Next let's look at the same method as above but for DPOAE's. We'll
%do this with a repeated measures ANOVA, with the effect of frequency and
%of time (but here only using two time points: baseline and 7+ post
%exposure time point).

%First set the names for the mice
clear; clc;
load('dataABR')
imagingAnimalsNE={'AEH147';'AEH155';'AEH162';'AEH165';'AEH168';'AEH170';'AEH174';'AEH178'};
HLexploreAnimals={'AEH141';'AEH144';'AEH145';'AEH149';'AEH150'};

%Now set the values
frequencyValues=[4 5.7 8 11.3 16 22.6 32 45.3 64]; 
imgDays={'day_pre';'day_1';'day_3';'day_5';'day_7';'day_10'}; 
nDays=size(imgDays,1); animals = [HLexploreAnimals;imagingAnimalsNE];
nFreq=length(frequencyValues); nDays=length(imgDays); nAnimals=length(animals);

%Store in a 3D matrix that is animals x frequency x day
thrDPOAE=nan(nAnimals,nFreq,nDays);

%Create the matrix
for ani=1:nAnimals
    animalStr=animals{ani,:};
    for dd=1:nDays
        dayStr=imgDays{dd,:};
        if ~isfield(dataABR.(animalStr),dayStr)
            continue
        end
        
        if isfield(dataABR.(animalStr).(dayStr),'both') && isfield(dataABR.(animalStr).(dayStr).both,'DPOAEthr')
            thrDPOAE(ani,1:nFreq,dd)=dataABR.(animalStr).(dayStr).both.DPOAEthr;
        end
    end
end

%Now we're going to select two specific time points, and only use
%frequencies that have been tested across all mice.

for i = 1:size(thrDPOAE,1) %Loop through mice
    nums1 = thrDPOAE(i,1:9,1); %Get baseline measurements
    nums2 = nan(1,9); day = size(thrDPOAE,3);
    while any(isnan(nums2(3:8)))
        nums2 = thrDPOAE(i,1:9,day);
        day = day - 1;
    end
    allValues(i,:) = [nums1 nums2];
end

%Now we want to format this into a table with the correct labels
for i = 1:18 %9 frequencies, 2 time points
    v = strcat('V',num2str(i));
    varNames{i,1} = v;
end

dataTable = array2table(allValues, 'VariableNames', varNames);

%Finally create a table that reflects the within subject factors, which
%here are frequency and time. 
factorNames = {'Frequency'; 'Day'};

freqLabels = repmat(frequencyValues(1:9),1,2)'; 
freqLabels = arrayfun(@num2str, freqLabels, 'UniformOutput', 0);

dayLabels = [ones(1,9) ones(1,9).*2]';
dayLabels = arrayfun(@num2str, dayLabels, 'UniformOutput', 0);

withinTable = table(freqLabels,dayLabels,'VariableNames',factorNames);

%Now that we have the tables set up, run the anova model
rmABR = fitrm(dataTable, 'V1-V18~1','WithinDesign',withinTable);
[rAnovaResults] = ranova(rmABR, 'WithinModel','Frequency*Day');
rAnovaResults

%Also do post-hoc tests for each frequency. Do as a paired t test, times
%two to correct for multiple comparisons.
for i = 1:9 %For 9 frequencies
    nums1 = allValues(:,i); nums2 = allValues(:,(i+9));
    [h,p] = ttest(nums1,nums2);
    disp(freqLabels{i}); disp(p*2); disp('');
end

%%
%Next is the statistics for 1D,1E (raw data also below)

%Fiirst is statistics on synapse counts. Below are the counts/IHC at each
%frequency measured.

clear; clc;

noise = [11.13 17.35 14.20 12.96 7.51 3.78;...
         12.22 16.35 18.51 10.23 8.67 6.95;...
         11.18 16.71 17.67 9.86 7.25 6.27;...
         11.93 15.97 19.41 8.08 5.26 3.42;...
         13.03 17.76 17.68 11.51 7.05 6.06;...
         12.20 16.75 18.09 10.17 6.63 5.75;...
         16.26 15.87 18.87 6.73 10.37 7.06];
     
sham = [11.23 16.55 19.68 16.18 12.80 14.76;...
        12.20 15.54 17.50 18.18 17.95 14.11;...
        8.02 14.45 19.92 18.48 16.66 11.82;...
        12.48 16.51 15.60 19.95 19.00 14.17];
    
allValues = [noise; sham];
    
%Now we want to format this into a table with the correct labels
for i = 1:6
    v = strcat('V',num2str(i));
    varNames{i,1} = v;
end

dataTable = array2table(allValues, 'VariableNames', varNames);

%Add in the exposure group variable
for i = 1:11
    if ismember(i,[1:7])
        expGroup{i,1} = 'Noise';
    else
        expGroup{i,1} = 'Sham';
    end
end
dataTable.expGroup = expGroup;

%Finally create a table that reflects the within subject factors, which
%here are frequency and time. 
factorNames = {'Frequency'};

freqLabels = [1:6]';
freqLabels = arrayfun(@num2str, freqLabels, 'UniformOutput', 0);

withinTable = table(freqLabels,'VariableNames',factorNames);

%Now that we have the tables set up, run the anova model
rmIHC = fitrm(dataTable, 'V1-V6~expGroup','WithinDesign',withinTable);
[rAnovaResults] = ranova(rmIHC, 'WithinModel','Frequency');
rAnovaResults

%Also do post-hoc tests for each frequency. Do as a two sample t test
%comparing noise to sham. One-sided test, Bonferonni corrected, so net no
%effect.

for i = 1:6 %For 6 frequencies
    nums1 = allValues(1:7,i); nums2 = allValues(8:end,i);
    [h,p] = ttest2(nums1,nums2);
    disp(freqLabels{i}); disp(p); disp('');
end

%%%%%%%%%%%%%%%
%Similar as above but now for the OHC counts, raw data is given below.

noise = [9.33 10.83 10.83 10.67 10 6.67;...
         10.17 10.50 10.50 10.17 9.50 7.50;...
         9.67 10.17 10.58 10.25 9.25 6.17;...
         10.00 10.33 10.00 10.33 9.33 3.67;...
         10.33 10.00 10.00 11.00 8.00 2.00;...
         10.00 11.00 10.83 11.00 10.33 8.00;...
         10.33 10.00 11.00 10.33 11.00 3.67];
     
sham = [10.00 10.67 10.00 10.67 9.67 9.67;...
        9.67 10.67 10.33 11.00 10.33 9.67;...
        9.67 10.33 10.67 11.00 11.00 9.67;...
        8.67 11.00 10.00 10.67 11.00 10.00];
    
allValues = [noise; sham];
    
%Now we want to format this into a table with the correct labels
for i = 1:6
    v = strcat('V',num2str(i));
    varNames{i,1} = v;
end

dataTable = array2table(allValues, 'VariableNames', varNames);

%Add in the exposure group variable
for i = 1:11
    if ismember(i,[1:7])
        expGroup{i,1} = 'Noise';
    else
        expGroup{i,1} = 'Sham';
    end
end
dataTable.expGroup = expGroup;

%Finally create a table that reflects the within subject factors, which
%here are frequency and time. 
factorNames = {'Frequency'};

freqLabels = [1:6]';
freqLabels = arrayfun(@num2str, freqLabels, 'UniformOutput', 0);

withinTable = table(freqLabels,'VariableNames',factorNames);

%Now that we have the tables set up, run the anova model
rmOHC = fitrm(dataTable, 'V1-V6~expGroup','WithinDesign',withinTable);
[rAnovaResults] = ranova(rmOHC, 'WithinModel','Frequency');
rAnovaResults

for i = 1:6 %For 6 frequencies
    nums1 = allValues(1:7,i); nums2 = allValues(8:end,i);
    [h,p] = ttest2(nums1,nums2);
    disp(freqLabels{i}); disp(p); disp('');
end

%%
%Next is the plotting for Figure 1I
load('example_run')

%Find all frequencies tested
for i = 1:length(Data)
    nums1(i) = Data{i}.cue.Speaker_1.Waveform.Frequency_kHz;
end
freqs = unique(nums1);

%Loop through trials in Data
allTrials = {};
for i = 1:length(freqs)
    allTrials{i} = NaN(length(Data),3);
end

for i = 1:length(Data)
    
    %Find the current intensity level and freqency
    ftype = Data{i}.cue.Speaker_1.Waveform.Frequency_kHz;
    ftype = find(freqs == ftype);
    allTrials{ftype}(i,1) = Data{i}.cue.Speaker_1.Level.dB_SPL;
    
    %Trial type
    if strcmp(Data{i}.Type,'adapt')
        allTrials{ftype}(i,2) = 1; %Regular trial
    else
        allTrials{ftype}(i,2) = 0; %Catch trial
    end
    
    %Define the outcome of the trial
    outcome = Data{i}.Result;
    if strcmp(outcome,'Go')
        allTrials{ftype}(i,3) = 1;
    else
        allTrials{ftype}(i,3) = 0;
    end
    
    
end

%Now eliminate all NaNs
for i = 1:length(allTrials)
    allTrials{i}(any(isnan(allTrials{i}), 2), :) = [];
end


%We will next plot the results
 
%Get the trials from the current frequency type
nums = allTrials{2};
for i = 1:length(nums)
    nums(i,4) = i;
end

%Decide the y axis limits
yMin = min(nums(nums(:,1)~=-1000,1)) - 5;
yMax = max(nums(nums(:,1)~=-1000,1)) + 5;

figure

%Plot regular trials
tmp = nums(nums(:,2) == 1,:);
tmp1 = tmp(tmp(:,3) == 1,:);
p1 = plot(tmp(:,4),tmp(:,1),'-ko','MarkerSize',9,'LineWidth',1);
hold on
p2 = plot(tmp1(:,4),tmp1(:,1),'ko','MarkerFaceColor','k','MarkerSize',9);

%Plot catch trials not -1000
tmp = nums(nums(:,2) == 0,:);
tmp = tmp(tmp(:,1) ~= -1000,:);
tmp1 = tmp(tmp(:,3) == 1,:);
p3 = plot(tmp(:,4),tmp(:,1),'ko','MarkerSize',9);
p4 = plot(tmp1(:,4),tmp1(:,1),'ko','MarkerFaceColor','k','MarkerSize',9);

%Plot catch trials that are -1000
tmp = nums(nums(:,2) == 0,:);
tmp = tmp(tmp(:,1) == -1000,:);
tmp1 = tmp(tmp(:,3) == 1,:);
p5 = plot(tmp(:,4),zeros(1,length(tmp(:,4)))+yMin,'o','MarkerSize',9,'Color',[0.5 0.5 0.5]);
p6 = plot(tmp1(:,4),zeros(1,length(tmp1(:,4)))+yMin,'o','MarkerFaceColor',[0.5 0.5 0.5],'MarkerSize',9,'Color',[0.5 0.5 0.5]);
ylim([yMin yMax])
xlabel('Trial Number','FontSize',16)
ylabel('Tone Level (dB SPL)', 'FontSize', 16)
title('8 kHz Track', 'FontSize',16)

legend([p1 p2 p3 p4 p5 p6],{'CS+ miss' 'CS+ hit' 'Range miss' 'Range hit' 'CS- withold' 'CS- FA'},...
         'FontSize',12)
hold off
box off

%%
%Next the plotting for Figure 1J

mouse = 12; %Since this is the 12th mouse in the structure
nTrialCutoff = 2; %Get rid of nonsense points
frequency = 1; %1 is 8 kHz, 2 is 32 kHz

%Load the data
load('exampleSessionData.mat')

%First plot the baseline data
i=4;
figure
for k = 1:2 %For two frequencies

if k == 1
    color1 = [0.5 0.5 1];
    color2 = [0 0 1];
else
    color1 = [1 0.5 0.5];
    color2 = [1 0 0];
end    

%Get the data
nTrials = Base.numTrial{mouse}{i}{k};
hits = nTrials.*Base.func{mouse}{i}{k};
thresh = Base.thresh{mouse}(i,k);

%Enforce the trial cutoff
index = nTrials > nTrialCutoff;
nTrials = nTrials(index);
hits = hits(index);
intensities = Base.trialType{mouse}{i}(index);


%Now do a logistic fit
xFit = [intensities(2)-5:0.25:intensities(end)+5];
[logitCoef,~] = glmfit(intensities(2:end),[hits(2:end)' nTrials(2:end)'],'binomial','logit');
yFit = glmval(logitCoef,xFit,'logit'); 

%Plot the results
plot(intensities(2:end),  hits(2:end)./nTrials(2:end),'Marker','.','LineStyle','none','MarkerSize',20,'Color',color1)
hold on
plot(xFit, yFit,'Color', color2, 'LineWidth', 1.5)   
plot([thresh thresh], [0 1],'Color',color2,'LineStyle','--')
title('Day -1')
xlim([0 100])
ylim([0 1])
set(gcf, 'Position',  [100, 100, 750, 500])
box off
set(gca, 'XTick', [0 25 50 75 100])
set(gca, 'YTick', [0 0.5 1])

%print out the FA rates
disp(num2str(hits(1)/nTrials(1)))

end
disp('-----')

%Now for the post exposure sessions
days = {'Day 0' 'Day 20'}; count = 1;
for i = [1 15] %Loop through the sessions
    
    figure
    
    for k = 1:2 %Loop through frequencies
        
    if k == 1
        color1 = [0.5 0.5 1];
        color2 = [0 0 1];
    else
        color1 = [1 0.5 0.5];
        color2 = [1 0 0];
    end   
    
    %Get the data
    nTrials = Expose.numTrial{mouse}{i}{k};
    hits = nTrials.*Expose.func{mouse}{i}{k};
    thresh = Expose.thresh{mouse}(i,k);
    
    %Enforce the trial cutoff
    index = nTrials > nTrialCutoff;
    nTrials = nTrials(index);
    hits = hits(index);
    intensities = Expose.trialType{mouse}{i}(index);
    
    
    %Now do a logistic fit
    xFit = [intensities(2)-5:0.25:intensities(end)+5];
    [logitCoef,~] = glmfit(intensities(2:end),[hits(2:end)' nTrials(2:end)'],'binomial','logit');
    yFit = glmval(logitCoef,xFit,'logit'); 
    
    %Plot the results
    plot(intensities(2:end),  hits(2:end)./nTrials(2:end),'Marker','.','LineStyle','none','MarkerSize',20,'Color',color1)
    hold on
    plot(xFit, yFit,'Color', color2, 'LineWidth', 1.5)  
    plot([thresh thresh], [0 1],'Color',color2,'LineStyle','--')
    title(days{count})
    xlim([0 100])
    ylim([0 1])
    set(gcf, 'Position',  [100, 100, 750, 500])
    box off
    set(gca, 'XTick', [0 25 50 75 100])
    set(gca, 'YTick', [0 0.5 1])
    
    
    %print out the FA rates
    disp(num2str(hits(1)/nTrials(1)))
    end
    disp('-----')
    count = count + 1;
end

%%
%Next is plotting for Figure 1K

%Get the threshold estimates
load('detectionBehaviorData.mat')


baseThresh = {[] []};
%Now align the baseline data correctly
for i = 1:length(Base.thresh) %Loop through mice
    threshs = Base.thresh{i};
    if i < 4
        threshs = [NaN NaN; threshs]; %Since these mice have only 3 baseline days
    end
    baseThresh{1} = [baseThresh{1} threshs(:,1)]; %8 kHz threshold data
    baseThresh{2} = [baseThresh{2} threshs(:,2)]; %8 kHz threshold data
end


expThresh = {[] []};
%Now align the exposure data correctly
for i = 1:length(Expose.thresh) %Loop through mice
    threshs = Expose.thresh{i};
    if i < 4
        threshs = [threshs(1:9,:);[NaN NaN];threshs(10:end,:)]; %These mice are missing Day 10
    end
    expThresh{1} = [expThresh{1} threshs(:,1)]; %8 kHz threshold data
    expThresh{2} = [expThresh{2} threshs(:,2)]; %8 kHz threshold data
end

%Now we can split the data into two groups

%Put baseline and post exposure together
baseThresh{1} = [baseThresh{1}; expThresh{1}];
baseThresh{2} = [baseThresh{2}; expThresh{2}];

%First define the 2 groups
sham = [5,6,7,8,11,13];
exp = [1,2,3,4,9,10,12];

means8sham = nanmean(baseThresh{1}(:,sham),2);
means32sham = nanmean(baseThresh{2}(:,sham),2);
er8sham = nanstd(baseThresh{1}(:,sham),[],2)./sqrt(length(sham));
er32sham = nanstd(baseThresh{2}(:,sham),[],2)./sqrt(length(sham));

means8exp = nanmean(baseThresh{1}(:,exp),2);
means32exp = nanmean(baseThresh{2}(:,exp),2);
er8exp = nanstd(baseThresh{1}(:,exp),[],2)./sqrt(length(exp));
er32exp = nanstd(baseThresh{2}(:,exp),[],2)./sqrt(length(exp));

figure
shadedErrorBar([1:13,15:2:25],means8sham,er8sham,'lineProps','-c');
hold on
shadedErrorBar([1:13,15:2:25],means32sham,er32sham,'lineProps','-m');
shadedErrorBar([1:13,15:2:25],means8exp,er8exp,'lineProps','-b');
shadedErrorBar([1:13,15:2:25],means32exp,er32exp,'lineProps','-r');
p1 = plot([1:13,15:2:25],means8sham,'-oc','lineWidth',2);
p2 = plot([1:13,15:2:25],means32sham,'-om','lineWidth',2);
p3 = plot([1:13,15:2:25],means8exp,'-ob','lineWidth',2);
p4 = plot([1:13,15:2:25],means32exp,'-or','lineWidth',2);
ax = gca;
xlabel('Day Re:Exposure', 'FontSize', 16)
ylabel('Threshold (dB SPL)','FontSize', 16)
ax.XTickLabel={'-5','0','5','10','15','20'};
legend([p1 p2 p3 p4],{'8 kHz Sham', '32 kHz Sham','8 kHz Exposed','32 kHz Exposed'},'FontSize',14)

%%
%Statistics for Figure 1K

%Get the threshold estimates
clear; clc;
load('detectionBehaviorData.mat')

baseThresh = {[] []};
%Now align the baseline data correctly
for i = 1:length(Base.thresh) %Loop through mice
    threshs = Base.thresh{i};
    if i < 4
        threshs = [NaN NaN; threshs]; %Since these mice have only 3 baseline days
    end
    baseThresh{1} = [baseThresh{1} threshs(:,1)]; %8 kHz threshold data
    baseThresh{2} = [baseThresh{2} threshs(:,2)]; %8 kHz threshold data
end


expThresh = {[] []};
%Now align the exposure data correctly
for i = 1:length(Expose.thresh) %Loop through mice
    threshs = Expose.thresh{i};
    if i < 4
        threshs = [threshs(1:9,:);[NaN NaN];threshs(10:end,:)]; %These mice are missing Day 10
    end
    expThresh{1} = [expThresh{1} threshs(:,1)]; %8 kHz threshold data
    expThresh{2} = [expThresh{2} threshs(:,2)]; %8 kHz threshold data
end

%Put baseline and post exposure together
baseThresh{1} = [baseThresh{1}; expThresh{1}];
baseThresh{2} = [baseThresh{2}; expThresh{2}];

%Now format a final matrix of all of the data
dataTable = [baseThresh{1}' baseThresh{2}'];
inds = find(isnan(dataTable(1,:))); dataTable(:,inds) = []; %Eliminate NaN's

%Now we want to format this into a table with the correct labels
for i = 1:length(dataTable)
    v = strcat('V',num2str(i));
    varNames{i,1} = v;
end
dataTable = array2table(dataTable, 'VariableNames', varNames);

%Add in the exposure group variable
for i = 1:13
    if ismember(i,[1,2,3,4,9,10,12])
        expGroup{i,1} = 'Noise';
    else
        expGroup{i,1} = 'Sham';
    end
end
dataTable.expGroup = expGroup;

%Next format the table containing the within subjects factors (the repeated
%measures).
factorNames = {'Frequency'; 'Day'};

freqLabels = [ones(17,1).*8; ones(17,1).*32];
freqLabels = arrayfun(@num2str, freqLabels, 'UniformOutput', 0);

dayLabels = vertcat([1:17]',[1:17]');
dayLabels = arrayfun(@num2str, dayLabels, 'UniformOutput', 0);

withinTable = table(freqLabels,dayLabels,'VariableNames',factorNames);

%Now that we have the tables set up, run the anova model
rmThresh = fitrm(dataTable, 'V1-V34~expGroup','WithinDesign',withinTable);
[rAnovaResults] = ranova(rmThresh, 'WithinModel','Frequency*Day');
rAnovaResults